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Abstract 

The advent of inexpensive, high-performance computers and new efficient algorithms 
have made possible the automatic recognition of numerically computed constants. In other 
words, techniques now exist for determining, within certain limits, whether a computed real 
or complex number can be written as a simple expression involving the classical constants 
of mathematics. 

These techniques will be illustrated by discussing the recognition of Euler sum con- 
stants, and also the discovery of new formulas for 7 and other constants, formulas that 
permit individual digits to be extracted from their expansions. 
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1. Introduction 

The advent of inexpensive, high-performance computers and new efficient algorithms 
have made possible the automatic recognition of numerically computed constants. In other 
words, techniques now exist for determining, within certain limits, whether a computed real 
or complex number can be written as a simple expression involving the classical constants 
of mathematics. 

The fundamental technique involved is that of finding integer relations. Let x = 


(%1,%2,°++,2,) be a vector of real numbers. z is said to possess an integer relation if 


there exist integers a; not all zero such that ajax, + dg%q +--+ + Gn%, = 0. By an inte- 
ger relation algorithm, we mean an algorithm that is guaranteed (provided the computer 
implementation has sufficient numeric precision) to recover the vector of integers aj, if it 
exists, or to produce bounds within which no integer relation can exist. 

The problem of finding integer relations among a set of real numbers was first studied 
by Euclid, who gave an iterative algorithm (the Euclidean algorithm), which when applied 
to two real numbers, either terminates, yielding an exact relation, or produces an infinite 
sequence of approximate relations. The generalization of this problem for n > 2 has been 
attempted by Euler, Jacobi, Poincare, Minkowski, Perron, Brun, Bernstein, among others. 
However, none of their algorithms has been proven to work for n > 3, and numerous 
counterexamples have been found. 

The first integer relation algorithm with the desired properties mentioned above was 
discovered by Ferguson and Forcade in 1977 [14]. In the intervening years a number of 
other integer relation algorithms have been discovered, including the “LLL” algorithm [16], 


the “HJLS” algorithm [15], and the “PSOS” [6] algorithm. 


2. The PSLQ Integer Relation Algorithm 

In 1991 a new algorithm, known as “PSLQ” algorithm, was developed by Ferguson 
[12]. It appears to combine some of the best features separately possessed by previous 
algorithms, including fast run times, numerical stability, numerical efficiency (i.e. suc- 


cessfully recovering a relation when the input is known to only limited precision), and a 


guaranteed completion in a polynomially bounded number of iterations. More recently a 
much simpler formulation of this algorithm was developed, and it has been extended to the 
complex number field [13]. This newer, simpler version of PSLQ can be stated as follows: 

Let x be the n-long input real vector, and let nint denote the nearest integer function 
(for exact half-integer values, define nint to be the integer with greater absolute value). 


Let y := \/4/3. Then perform the following: 
Initialize: 
1. Set the n x n matrices A and B to the identity. 


2. For k := 1 to n: compute sz := fe ore endfor. Set t = 1/s,. For k := 1 to n: 


YR i= txy; Sp := ts,; endfor. 


3. Compute the n x (n — 1) matrix A as follows: 
For i:= 1 ton: forj :=i+1ton—1: set H,; := 0; endfor; if i <n—1 then set 
Ay, = 8441/8;; for J = 1 toi — 1: set Ay; := —yy;/(8;8j;41); endfor; endfor. 

4. Perform full reduction on H, simultaneously updating y, A and B: 
Fort = 2 ton: for 7 <=4—ltol step —1>t <= nint( Ay / 7); 2 = uf is for k= 1 
10-4: Ay = Bet gs endior: tor k=) ton: Ay Ag —tAje, Big = Bap +t Be; 


endfor; endfor; endfor. 
Repeat until precision is exhausted or a relation has been detected: 
1. Select m such that 7'|H;;| is maximal when i = m. 


2. Exchange entries m and m+1 of y, corresponding rows of A and H, and corresponding 


columns of B. 


3. Ifm <n-—2 then update H as follows: 


Set to := He + is ee ty := Anm/to and to := Hmmii/to. Then for 7 := m to 
nis =a te = EB ae ae a ti = Stole + ita enidior: 
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4. Perform block reduction on H, simultaneously updating y, A and B: 


For 1:= m+1 ton: for j := min(t—1,m+1) to 1 step —1: ¢ := nint(H,;/H;;); y= 
y; tty; for k := 1 to 7: Hy := Hy —tHjx; ; endfor; fork := lton: Ay := 
Aix — tAjr, Bez = Bry +tByi; endfor; endfor; endfor. 


5. Norm bound: Compute M := 1/max;|H;|, where H; denotes the j-th row of H. 


Then there can exist no relation vector whose Euclidean norm is less than M. 


6. Termination test: If the largest entry of A exceeds the level of numeric precision 
used, then precision is exhausted. If the smallest entry of the y vector is less than 
the detection threshold, a relation has been detected and is given in the corresponding 


column of B. 


With regards to the termination criteria in step 6, it sometimes happens that a relation 
is missed at the point of potential detection because the y entry is not quite as small as 
the detection threshold being used (the threshold is typically set to the “epsilon” of the 
precision level). When this happens, however, one will note that the ratio of the smallest 
and largest y vector entries is suddenly very small, provided sufficient numeric precision is 
being used. In a normal computer run using the PSLQ algorithm, prior to the detection 
of a relation, this ratio is seldom smaller than 10~?. Thus if this ratio suddenly decreases 
to a very small value, such as 10~?°, then almost certainly a relation has been detected 
— one need only adjust the detection threshold for the algorithm to terminate properly 
and output the relation. When detection does occur, this ratio may be thought of as a 
“confidence level” of the detection. 

As a general rule, one can expect to detect a relation of degree n, with coefficients 
of size 10”, provided that the input vector is known to somewhat greater than mn digit 
precision, and provided that computations are performed using at least this level of numeric 


precision. 


3. Applications of the PSLQ Algorithm 

There are a number of applications of integer relation detection algorithms in compu- 
tational mathematics. One application is to analyze whether or not a given constant a, 
whose value can be computed to high precision, is algebraic of some degree n or less. This 
can be done by first computing the vector x = (1,a,a,---,a@”) to high precision and then 
applying an integer relation algorithm to the vector x. If a relation is found, this integer 
vector is precisely the set of coefficients of a polynomial satisfied by a. Even if a relation is 
not found, the resulting bound means that a@ cannot possibly be the root of a polynomial 
of degree n, with coefficients of size less than the established bound. Even negative results 
of this sort are often of interest. 

We have performed several computations of this type [6]. These computations have 
established, for example, that if Euler’s constant y satisfies an integer polynomial of degree 
50 or less, then the Euclidean norm of the coefficients must exceed 7 x 10!”. Computations 
of this sort have also been applied to study a certain conjecture regarding the Riemann 


zeta function. It is well known [10] that 


ee 1 
(2) = 3)0 
a 2 es 
5 (oe) (-1 k-1 
SC ame > 
2% BE) 
BG ah 
a Say 
1 a ee, 
These results have led some to suggest that 


Fi 0 Say 


might also be a simple rational or algebraic number. Unfortunately, integer relation calcu- 


lations [3] have established that if Z; satisfies a polynomial of degree 25 or less, then the 


Euclidean norm of the coefficients must exceed 2 x 10°”. 


4. Euler Sums 


In response to a letter from Goldbach, Euler considered sums of the form 


(1 | = bass oa) (eae 


k=1 


Euler was able to give explicit values for certain of these sums in terms of the Riemann 

zeta function. For example, Euler found an explicit formula for the case m = 1, n > 2. 

Little progress has been made on this problem in the intervening years, although special 

cases of Euler’s results have been rediscovered numerous times (see [7] for some references). 
In April 1993, Enrico Au-Yeung, an undergraduate at the University of Waterloo, 

brought to the attention of Prof. Jonathan Borwein the curious fact that 

3 (1+ ; foes =): k-? = 459987... ea) = le 


k=1 


based on a computation to 500,000 terms. Borwein’s reaction was to compute the value of 
this constant to a higher level of precision in order to dispel this conjecture. Surprisingly, 
a computation to 30 and later to 100 decimal digits still affirmed it. 

Intrigued by this empirical result, numerical values were computed for several of these 
and similar sums, which were termed Euler sums. The resulting values were analyzed 
using integer relation searches. These efforts produced even more empirical evaluations, 
suggesting broad patterns and general conjectures. Ultimately proofs were found for many 
of these experimental results. 

We will consider here only the following two classes of Euler sums. Other classes are 


studied in [4]. 


(oe) 1 m 

som): = ae ; ) (k+1)™ nm > 1,1 > 2, 
te 1 (=<1)*4 m 

8a(m, Nn) = » (1-5 a k (k+1)™ mel,n>2, 
k=1 


Explicit evaluations of some of the constants in these classes are presented with proofs in 


[8] and [9]. 


5. Numerical Techniques 

It is not easy to compute numerical values of any of these Euler sums to high preci- 
sion. Straightforward evaluation using the defining formulas, to some upper limit feasible 
on present-day computers, yields only about eight digits accuracy. Because integer rela- 
tion algorithms require much higher precision to obtain reliable results, more advanced 
techniques must be employed. 

Our approach to computing numerical values of these sums involves the compound 
application of the Euler-Maclaurin summation formula (see [2, p. 289]), which can be 
stated as follows. Suppose f(t) has at least 2p + 2 continuous derivatives on (a,b). Let 
D be the differentiation operator, let B, denote the k-th Bernoulli number, and let B,(-) 


denote the k-th Bernoulli polynomial. Then 


Sia) = f Float + SU@) + FO) 


SID F(0) — DY Fla] + Fylasd) () 


where the remainder R,(a, b) is given [2, p. 289] by 


Ryle!) = Gay ff Bonsall — EDP) ae 


(2p + 2 
We will briefly present a method for computing s;,(m, 7). See [4] for more details. Let 
h(k) = O4_, 1/7 and f(t) = 1/t. By applying the Euler-Maclaurin summation formula, it 


can be seen that 


1 1 1 1 1 
hk) = y+ Ink + 5 — aoe + Googe ~ aeons + aa0R8 
i 1 1 1 
ee — pues i + O(k7'8). (2) 


~ 132k10 " 32760k!2 2h 816016 


We will use h(k) to denote this particular approximation (i.e., (2) without the error 


term). Now consider the sum 


Let c be a large integer, and let g(t) = h™(t)(t+1)-". Applying the Euler-Maclaurin 


summation formula (1) again, we can write 


sn(m,n) = (1 ; pe =) a ae 


es (k+1) 
= Hames +f goat + 5ole+) 
os age e+) + O(c). (3) 


This formula suggests the following computational scheme. First, explicitly evaluate 
the sum 7¢_, h™(k)(k +1) for c = 10%, using a numeric working precision of 150 digits. 
Secondly, perform the symbolic integration and differentiation steps indicated in formula 
(3). Finally, evaluate the resulting expression, again using a working precision of 150 digits. 
The final result should be equal to s,(m,n) to approximately 135 significant digits. 

We have performed many computations of this type. The integration and differentiation 
operations required in (3) can be handled using a symbolic mathematics package, such as 
Maple [11] or Mathematica [17]. The explicit summation of the first c terms, as indicated 
in (3), could be performed by utilizing the multiple precision facility in the Maple or 
Mathematica packages. However, it was found that the MPFUN multiple precision package 
and translator developed by one of us [3] was significantly faster for this purpose. 

Whatever software is used, this explicit summation is an expensive operation. For 
example, the evaluation of s,(3,4) to 10° terms, using the MPFUN package with 150- 
digit precision arithmetic, requires 20 hours on a “Crimson” workstation manufactured by 
Silicon Graphics, Inc. Thus while such runs can be made, clearly this is pressing the limits 
of current workstation technology. Fortunately, it is possible to perform such computations 
on a highly parallel computer system. The details of this parallel algorithm are given in 


[4]. 


6. Application of PSLQ to Euler Sums 
The problem of recognizing Euler sum constants is well suited to analysis with integer 


relation algorithms. We will present but one example of these computations. Consider 


co 1 (—1)**1 2 7 

sa(2,3) = S (1-544 , (k+1)3 
k=1 

0.156166933381176915881035909687988193685776709840 - - - 


Based on experience with other constants, we conjectured that this constant satisfies a 
relation involving homogeneous combinations of ¢(2),¢(3), ¢(4), ¢(5), In(2), Lig(1/2) and 
Lis(1/2), where Li,,(2) = 7%, 2*k-” denotes the polylogarithm function. 

The set of terms involving these constants with degree five (see section 7) are as follows: 
Lis(1/2), Lia(1/2) n(2), In°(2), ¢(5), ¢(4)In(2), (3) n2(2), ¢(2)In*(2), ¢(2)¢(3). When 
Sa(2,3) is augmented with this set of terms, all computed to 135 decimal digits accuracy, 
and the resulting 9-long vector is input to the PSLQ algorithm, it detects the relation 
(480, —1920, 0, 16, 255, 660, —840, —160,360) at iteration 390. Solving this relation for 


8a(2,3), we obtain the formula 


s4(2,8) = ALis(1/2) — 3 In°(2) ~ 32¢(5) — C(4) n(2) + £6(8) m2(2) 
+36(2) n°) — $6(2)¢(3) 
= 4Lis(1/2) — 5 n*(2) — =6(6) - aot In(2) + (3) In?(2) 
2 3 3 2 
Hage In°(2) — vid ¢(3) 


(recall that ¢(2n) = (27)?"|Bon|/ [2(2n)!]). 
When the relation is detected, the minimum and maximum y vector entries are 1.60 x 
10-184 and 5.98 x 10~?°, respectively. Thus the confidence level of this detection is on the 


order of 107! 


, indicating a very reliable detection. 
Although 135-digit input values and 150-digit working precision were used when this 
relation was originally detected, the fact that the maximum y-vector entry is only 10779 at 


detection implies that such high levels of numeric precision are not required in this case. 


Indeed, the above relation can be successfully detected using only the 50-digit input values 
listed above and 50-digit working precision when performing the PSLQ algorithm. 

Table 1 lists a number of the formulas that have been found by this procedure. Table 
2 lists some relations between the s, and s, constants. Others of both classes be found in 
[4]. It should be emphasized that the results in Tables 1 and 2 are not established in any 
rigorous mathematical sense by these calculations. However, in each case the “confidence 
level” (see section 3) of these detections is less than 10~°°, and in most cases is in the 


neighborhood of 10710, 


7. New Formulas for 7 and Related Constants 

Many readers may be already familiar with the recent paper by the authors and Peter 
Borwein that gives a new method for computing individual digits of 7 and certain other 
related constants. A brief review of these results is as follows. First, we inquire whether 7 


satisfies a relation of the form 


_ a ao | ay, ag a3 a4 
rs > i6i Ee Skil Sk+2° 843 8k 
| Os | 46 | CE 
'8kK+5  8k+6  8k+7 


where a; are rational numbers. Indeed it does. Such a formula can be found by separating 
the right hand side of the above expression into eight summations, numerically evaluating 
each to high precision, appending the numerical value of 7, and applying PSLQ to the 


resulting 9-long vector. When this is done, PSLQ discovers the following formula: 


mo 17 4 2 1 1 
= - = - 4 
. » 16 Fes 8k+4 8k+5 8k+6 (4) 


A similar formula was discovered by Ferguson. These two formulas form the basis of a 
two-dimensional lattice of formulas of this form. 
The significance of these formulas for the computation of 7 can be seen as follows. Let 


S; be the first of the sums in the above formula (4). Then we can write 


oo 16¢-* 
frac(167S;) = rg) 


k=0 


(mod 1) (5) 
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su(3,2) = 2C(5) +¢(2)¢(3) 
sn(3,3) = ——26(6) + 26703) 
su(3,4) = Fec(1) — Pctaye(a) + 2¢(2)6(5) 
su(3,6) = SEC(9) — c(ayets) — Lc(3)¢(6) + C8) + 36(2)¢(7) 
su(4,3) = —Sc¢7) + Fe(aycta) — 5¢(2)¢(5) 
su(5.4) = Ec(9) + 686(4)¢(5) - SPc(3)6(6) — 5¢4(8) + c(2)«I7) 
sn(6,3) = —Ss-c(9) — 248¢(4)¢(5) + “c(aye(o) + FZeh(ay - Sccapeva) 
su(72) = SPE c(o) 4 PE ccayc(sy + = c(a)e(6) + 56¢%(8) 
+S Zeayciz) 
8a(2,2) = 6Li4(1/2) + + m'(2) — = 64) + 5¢(2) In*(2) 
s4(2,8) = ALig(1/2) — 3 n°(2) — 3£¢(5) ~ “26(4) (2) + £¢(8) m2(2) 
+ 5¢(2)Im(2) — 36(26(8) 
8a(3,2) = —24Li5(1/2) + 6In(2)Li,(1/2) + = m2) + (5) - (4) In(2) 


5 ; f 
+ 56(2) n*(2) + 56(2)¢(3) 


Table 1: Experimentally Detected Formulas 


845495), (1,7) + 2114685), (2,6) + 148902s),(3, 5) — 13360s,(4, 4) — 19785;,(5, 3) 
~2718587s,,(1,8) — 1645256646), (2, 7) — 178042944s,(3, 6) — 88947862s),(4, 5) 
+3863940s),(5, 4) + 6721005), (6, 3) 

—142694085;,(1,9) + 2578470s),(2, 8) + 2815376s,(3, 7) + 58145505, (4, 6) 
+62388845),(5, 5) + 3938912s,(6, 4) + 1122784s,(7, 3) — 1860s,(8, 2) 
-63164285C (10) 

321¢(10) — 440¢°(5) — 720¢(3)¢(7) — 80¢7(3)¢(4) + 560¢(2)¢ (3)¢(5) 

—40s1,(2, 8) + 160s,(3, 7) 

—1691755503s;(1, 10) — 3172589688s),(2, 9) + 837511504), (3, 8) 


—7302717576s;,(4, 7) — 13958660016s,(5, 6) — 12910466064.) (6, 5) 
—7099332912s;,(7, 4) — 17732126885;,(8, 3) + 6583605) (9, 2) 
+53491434679¢(11) — 21868248971¢(2)¢(9) 

—589¢(11) + 322¢(5)¢(6) + 756¢(4)¢(7) + 254¢(3)¢(8) — 336¢2(3)¢(5) 
—368¢(2)¢(9) + 80¢(2)¢3(3) — 16s;,(3, 8) — 485,,(4, 7) 

1152s,(2, 4) + 640s, (3, 3) — 7680 In(2)Lis(1/2) + 64 In®°(2) — 1881¢(6) 
+7440¢(5) In(2) — 1680¢(4) In?(2) — 1120¢(3) In?(2) + 864¢(3)¢(3) 
—640¢ (2) In4(2) — 432¢(2)¢(3) In(2) 


Table 2: Experimentally Detected Relations 
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216%" (mod 8k41). » 2k. 6s * 
= | dl 6 
» 8k +1 2» ee oe) (6) 


k=0 k=d+1 


where frac denotes fractional part, i.e. the value mod 1. 

The numerator of the first summation in (6) can be rapidly evaluated by means of the 
binary algorithm for exponentiation, where each operation is performed modulo the integer 
8k+1. These calculations can be done with either integer or floating-point arithmetic, pro- 
vided the format being used has enough accuracy to exactly represent the integer d?. Once 
an individual exponentiation operation is complete, the resulting integer value is divided 
by 8k + 1, using floating-point arithmetic, and added to the sum modulo 1. Only a few 
terms are required of the second summation in (6), since they rapidly become smaller than 
the “epsilon” of the floating-point arithmetic system being used. The resulting fractional 
value, when expressed in base 16 notation, gives the hexadecimal digits of 7 beginning at 
position d+ 1. 

Here are a number of other formulas of this type. As before, these formulas were 


originally found using PSLQ searches. 


» 1a 1 144 216 72 B40 
a < 3 bak on — (6k +2)? (6k +3)? (6k +4)? ° (k+5)? 
; OY 16 16 8 16 
= dX 16 a +12 (k+2) (8k+3)2 (844)? 
4 Bit 8 
(8K +5)2 (8k +6)2  (8k+ =a 
1 FT 8 16 40 8 28 
log'(2) = 6 6k Lie ' e412 (@k+22 Bk+32 Bk+42 
4 i 2 
"(8k +5)? (8k+6)2 (8k+ aa 


Full details of these calculations, as well as formal proofs of the above formulas, can be 


found in [5]. 


8. A General Constant Recognition Procedure 
In all of the cases mentioned above, the authors of the respective studies had “hunches” 


beforehand as to what form the resulting formulas might take. Frankly, some insight of 
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sort is invaluable in avoiding what otherwise is a combinatorial explosion in the number 
of possible terms. It simply is not possible to perform integer relation searches with every 
conceivable term. In fact, if the constant is known to only limited precision, the number 
of terms that can be considered in an integer relation search may be limited to a handful. 

Nonetheless, it does appear feasible to define general search procedures that are often 
successful in recovering the analytic form of many constants that naturally appear in 


mathematical calculations. The authors present the following procedure as an example: 


1. Using PSLQ and full precision, check if a’ is algebraic of degree n, for j up to m. 


2. Using PSLQ and full precision, check if a is given by a multiplicative formula of the 


form 


0 = a,log(a) + ag log(2) + a3 log(3) + a4 log(5) + --- + a, log(p,.) 


+G4r41 log(c1) a Ar+2 log(c2) ae Po Ar+s log (cs) 


where p; is the k-th prime, and where cy are a selected set of transcendentals. 


3. Using PSLQ and quad precision, check if @ is given by a linear formula of the form 


0 = ao + a1 + A2g%1 + a3X%q+°--4 At41%t 


where 21, %2,-+-,2,; are each a product of up to three constants from a selected set 


of algebraic and transcendental constants. 


4. If a tentative relation is found in the previous step using quad precision, then check 


it using full precision. 


In the above, m,n,r,s,t are parameters that can be set to adjust the amount of computer 
time one is willing to expend in the search. The authors have tried, for example, m = 
Sw 32, ¥ = 20. 6 = 20; and f= 3. 


Some examples of constants recognized by above procedure are the following: 
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. The root near 1.38851367 of the polynomial 


—14250992566272 + 1934517374976t° — 37548447232t!? + 3072863296t'8 


+3789095144t74 — 408473063t°° — 43879700t°6 — 5815353¢*? 


S3TOGLi = 7Os8tee = Bolk Sad 190, 0a aa 


. The definite integral 


i (res dh = ane 
; 16r(1/4) 


. The definite integral 


m/4 ¢? dt 
= = 1 In(2) /4 
[ ag) 72/16 + rln(2)/4+G 


where G denotes Catalan’s constant 


. The definite integral 


2 


2_1)(t4+1) —-:16(2+. 2) 


1 ¢? In(t) dt 7 
he 
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